/*
 * l2lossl1reg.cpp
 *
 *  Created on: Feb 28, 2011
 *      Author: tqlong
 */

#include "perceptron.h"
#include <armadillo>
#include <algorithm>
#include <limits>
using namespace arma;

#define DOUBLE_INFINITY std::numeric_limits<double>::infinity()

void L1LossL2Reg::setParameters(const std::vector<double>& params)
{
  if (params.size() > 0) alpha_ = params[0];
  if (params.size() > 1) C_ = params[1];
  if (params.size() > 2) maxiter_ = (int)params[2];
}

/** Solving w \in Rd
 *      min_w 1/2 dot(w,w) + C \sum xi_i
 *      s.t. xi_i >= 0 && xi_i >= 1-y_i dot(w,x_i)
 *  in dual form a \in Rn
 *      min_a 1/2 dot(a,Qa) - dot(1,a)
 *      s.t. 0 <= a <= C
 *      where Q_ij = y_i y_j dot(x_i,x_j)
 */
double L1LossL2Reg::update(Dataset& data)
{
  a_ = zeros<vec>(data.n(), 1);
  w_ = zeros<vec>(data.dim(), 1);
  prepareQD(data);

  int iter = 0;
  double error = 1;
  std::cout << "maxiter = " << maxiter_ << "\n";
  while (iter < maxiter_) {
    iter++;

    vec a_old(a_);
    iterate(data);

    error = 0;
    for (int i = 0; i < data.n(); i++) {
      double pred = predict(data.x(i));
      error += (pred * data.y(i) <= 0 ? 1 : 0);
      //std::cout << "err " << i << " " << data.y(i) << " " << pred << " "<< (pred * data.y(i) <= 0) << "\n";
    }
    std::cout << "iter = " << iter << " error = " << error << "\n";

    double a_change = arma::norm(a_old-a_,"inf");
    std::cerr << "iter = " << iter << " error = " << error
        << " change = " << a_change << "\n";
    if (Compare(1e-5).equal(a_change, 0.0)) {
      std::cout << "Changes too small\n";
      break;
    }
  }
  return error;
}

void L1LossL2Reg::iterate(Dataset& data)
{
  Compare comp;
  for (int i = 0; i < data.n(); i++) {
    double G = data.y(i)*dot(w_, data.x(i)) - 1 + Dii_[i]*a_[i];

    double PG = G;
    if (comp.equal(a_[i], 0.0))
      PG = std::min(G, 0.0);
    else if (comp.equal(a_[i], U_))
      PG = std::max(G, 0.0);

    if (!comp.equal(PG, 0.0)) {
      double a_old = a_[i];
      a_[i] = std::min(std::max(a_[i] - G/Qii_[i], 0.0), U_);
      w_ += data.y(i) * (a_[i]-a_old) * data.x(i);
    }
  }
}

double L1LossL2Reg::predict(const arma::vec& x)
{
  double y = dot(w_, x);
  predict_ = y > 0 ? 1 : -1;
  return predict_;
}

void L1LossL2Reg::prepareQD(Dataset& data)
{
  std::cout << "prepare Q and D for L1Loss L2Reg\n";
  Qii_ = zeros<vec>(data.n(), 1);
  Dii_ = zeros<vec>(data.n(), 1);
  U_ = C_;
  for (int i = 0; i < data.n(); i++)
    Qii_[i] = dot(data.x(), data.x()) + Dii_[i];
}

void L2LossL2Reg::prepareQD(Dataset& data)
{
  std::cout << "prepare Q and D for L2Loss L2Reg\n";
  Qii_ = zeros<vec>(data.n(), 1);
  Dii_ = ones<vec>(data.n(), 1) * 0.5/C_;
  U_ = DOUBLE_INFINITY;
  for (int i = 0; i < data.n(); i++)
    Qii_[i] = dot(data.x(), data.x()) + Dii_[i];
}
